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ABSTRACT 

The Space-Station uses small rocket motors, called thrusters, for ori- 
entation control. Because of the lack of viable design tools for small 
rockets, the initial thruster design was basically a very small version of a 
large rocket motor. Thrust measurements of the initial design were lower then 
predicted. To improve predictions it was decided to develop a version of the 
RPLUS2D reacting flow code for thruster calculations. RPLUS2D employs an imp- 
licit finite volume, lower-upper symmetric successive overrelaxation (LU-SSOR) 
scheme for solving the complete two-dimensional Navier-Stokes equations and 
species transport equations in a coupled and very efficient manner. The com- 
bustion processes are modeled by a 9-species, 18 step finite-rate chemistry 
model, and the turbulence is simulated by a Baldwin-Lomax algebraic model. 

The code in this work is extended to handle multiple subsonic inlet conditions 
where the total mass flow is governed by conditions calculated at the 
thruster-throat. Results are shown for a thruster design where the overall 
mixture ratio is hydrogen rich. A calculation of a large area ratio divergent 
nozzle is also presented. 


INTRODUCTION 

For the space station relatively small rocket engines, called thrusters, 
are needed for attitude control. Using current successful designs for large 
rockets and the TDK code (which is based on the method of characteristics) a 
thruster was designed, built and tested. Thrust measurements of the small 
rocket motor were about 10 percent lower than predicted. As the thrusters may 
be expected to operate for long periods of time this error can have a large 
effect. In addition, it is very costly and time-consuming to design, fab- 
ricate, and test rocket engines. To cut down on costs more accurate predic- 
tions are needed. Several sources of error were identified. First, the TDK 
code does not predict the boundary layer as accurately as needed for thrusters. 
The TDK code calculates the flowfield in two passes. In the first pass, no 
boundary layer is calculated, then the flow is corrected for a boundary layer 
in the second pass. The boundary layer has a much larger effect in the rather 
small thruster, so calculating the boundary layer in this manner is not accu- 
rate enough. Secondly, the TDK code separates the flow into separate flow- 
streams. Unless a large number of flowstreams is specified, no mixing is 
allowed between flowstreams. This causes a large demand on computer resources. 


Finally, the TDK code uses an equilibrium chemistry model, which will tend to 
overpredict mixing and the extent of chemical reactions. Rather than adjust 
the constants in the TDK code to correctly predict the flow fields in 
thrusters, it was decided to develop a more efficient modern code. After a 
period of study, the RPLUS2D code was chosen for development. RPLUS2D has 
previously been developed for the study of mixing and chemical reactions in 
flow fields of ramjets and scramjets. 

RPLUS2D employs an implicit finite volume, lower-upper symmetric over- 
relaxation (LU-SSOR) scheme for solving the complete two-dimensional Navier- 
Stokes equations and species transport equations in a coupled manner, with 
real gas properties. The LU-SSOR technique was originally proposed by Jameson 
and Turkel (ref. 1) and developed by Yoon and Jameson (ref. 2). The code was 
extended to treat chemically reacting flows by Shuen and Yoon (ref. 3). The 
LU-SSOR scheme treats the inviscid flux and chemical source terms implicitly, 
but the viscous term explicitly. Despite being implicit, the LU scheme req- 
uires only scalar inversions for nonreacting flows and block-diagonal inversion 
for reacting flows, which saves considerable CPU time without using a diag- 
onal ization procedure. Furthermore, because of the special structure of the 
left-hand side (LHS) LU operators in the LU-SSOR scheme, the two L and U 
operators can be fully vectorized while the LHS operators in most other 
implicit schemes are not fully vectorizable. 


GOVERNING EQUATIONS 

The two-dimensional Navier Stokes and species transport equations for a 
chemically reacting flow with N s species, in Cartesian coordinates, are 
given by equations (1) and (2). 
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The diffusion velocities are found by Fick's law 
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where Di m is the binary diffusivity for species i in the gas mixture. 

The temperature and pressure are calculated iteratively from the following 
equations : 
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Here Yj, X\ and M^ are the mass fraction, molar fraction and 
molecular weight for species i. 


THERMODYNAMIC AND TRANSPORT MODELS 

The specific heat, thermal conductivity and viscosity for each species 
are determined by fourth-order polynomials of temperature, such as: 


Cp^ = A i + BjT + CjT 2 + DjT 3 + EjT 4 
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The coefficients of these polynomials are supplied by McBride (ref. 4) and 
are valid up to a temperature of 6000 K. The specific heat of the gas mixture 
is obtained by concentration weighting of each species. The thermal conduc- 
tivity and viscosity of the mixture, however, are calculated using Wilke's 
mixing rule (ref. 5) . 

The binary mass diffusivity Dj, between species i and j is 
obtained using the Chapman-Enskog theory in conduction with the Lennard-Jones 
i ntermol ecul ar potential-energy functions (ref. 5). The diffusion of a 
species in the gas mixture is approximated by Fick's Law, i.e., treating the 
species i and the surrounding gas as a binary gas mixture, and the diffusion 
velocity for each species is calculated using equation (2). 


CHEMISTRY MODELS 


In the present study, the finite-rate reaction of hydrogen and air are 
considered for the chemistry model. The reaction model of Brabbs (ref. 6) is 
used. This model employs 9 species and 18 elementary reactions. Nitrogen is 
treated as an inert species. Nitrogen was not specified for the thruster 
calculation. For a set of Nr elementary reactions involving N species, 
the rate equations can be written in the following general form: 
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Here Cj is the molar concentration of the jth species. The forward reaction 

rate constant kf for the i th elementary reaction is given by the 
i 

Arrhenius equation: 

Bi 

kf - AfT exp (-Ef/RT) (7) 

where E} represents the activation energy, and A} and Bj are con- 
stants. The backward reaction rate constant is calculated using the forward 
reaction rate constant and the equilibrium coefficient k^ , 
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The rate of change of molar concentration of species j by reaction 
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The total rate of change of molar concentration of species j is then 
given by 


• Nr . 
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the chemical source terms can be found from 


Wj - C jMj 


( 11 ) 


TURBULENCE MODEL 


The zero-equation algebraic turbulence model developed by Baldwin-Lomax 
(ref. 7) is used. This model is a two-layer model in which an eddy viscosity 
is calculated for an inner and outer region. In both of the regions the 
distribution of vorticity is used to determine the length scales. The primary 
advantage of this model is its simplicity and that boundary layer thickness 
does not need to be calculated. 


NUMERICAL TECHNIQUE 

Many numerical techniques have been used to solve the set of equations 
governing chemically reacting flows. Among these techniques, the explicit 
schemes are generally slow in convergence when the flow involves high rates of 
heat release or zones of recirculation. Most implicit schemes, on the other 
hand, require the inversion of block matrices and become exceedingly expensive 
when the chemical system involves a large number of species. In the present 
study, the lower-upper symmetric overrelaxation (LU-SSOR) scheme of Yoon and 
Jameson 2 is adopted to solve the two-dimensional Navier-Stokes and species 
transport equations. The LU-SSOR scheme employs implicit Newton iteration 
technique to solve the finite-volume approximation of the steady-state version 
of the governing equations. The convergence of the Newton iteration method is 
assured by the diagonal dominance of the coefficient matrices of the LU-SSOR 
scheme. The detailed derivations of the LU-SSOR scheme for nonreacting and 
reacting flows can be found in references 2 and 3, repsecti vely. The final 
discretized form of the equations for chemically reacting flows are given by 
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where the Jacoabian matrices 
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and T 

are defined by 
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Here because of the stiffness of the equations for reacting flows, the 
chemical source terms have been treated implicitly. Bussing and Murman 
(ref. 8), and Eklund et al . (ref. 9) have shown that implicit treatment of the 
chemical source terms is equivalent to rescaling the governing equations in 
time and thus stablizes the numerical scheme. 

In equation (12), D x ~ and Dy - are backward difference operators and D x + 

+ 

and Dy are forward difference operators. A + , A", B + and B~ are 
constructed so that the eigenvalues of "+" matrices are nonnegative and those 
of matrices are nonpositive, as follows: 


A + = 1/2 (A + 7 A I) 
A- = 1/2 (A - 7 A I) 
B+ = 1/2 (B + 7 B D 
B- - 1/2 (B - 7 B I> 


where 
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Here, A/\ ancl *B represent eignevalues of Jacobian matrices. 

After expanding out the difference operators on the LHS, equation (12) 
becomes 
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It is interesting to note that for nonreacting flows (S = 0 and T = 0 in 
eq. (16)) the present numerical method eliminates the need for banded block 
matrix inversions without using a diagonal ization procedure. In fact, with a 
forward diagonal sweep for the inversion of the first operator and a backward 
diagonal sweep for the inversion of the second operator, only scalar diagonal 
inversions are needed to solve equation (16) for nonreacting flow problems. 

For reacting flows, due to the presence of the chemical source Jacobian 
T, the first operator on the left-hand side of equation (16) now requires 
block diagonal inversions. However, since in the present formulation the flow 
equations (continuity, momentum, and energy equations) have no source terms. 


6 


the first four rows of the diagonal block of this operator have nonzero terms 
only in the diagonal. As a result, the first operator of equation (16) can be 
inverted in essentially two separate steps: the scalar diagonal matrix 

inversion for the flow equations and the block diagonal matrix inversion for 
the species equations. In contrast, most other implicit schemes require block 
banded (tridiagonal or pentadiagonal ) matrix inversions for both operators for 
either the reacting or nonreacting flows. Here, two advantages in efficiency 
are obtained by the LU-SSOR scheme over other implicit schemes: (a) the block 

size is smaller in the LU-SSOR scheme. Since the operational count is prop- 
ortional to for inverting a N x N matrix, the saving in CPU time can be 
significant; (b) the operational count for inverting a block diagonal matrix 
is only a fraction of that for inverting a block banded matrix. 

Another important efficiency advantage of the LU-SSOR scheme is the vec- 
torizability of its implicit operators. For example, the inversion of the two 
LHS operators of equations (16) can be fully vectorized if the solution sweeps 
are taken along the diagonal directions (the directions that i + j increases 
or decreases). It is noted that the vector i zabi 1 i ty is rather unique to the 
LU scheme, and, few other implicit scheme can have the LHS operators fully 
vectorized. The vectori zabi 1 i ty of the LHS operators is particularly sig- 
nificant for reacting flow calculations, since the most time-consuming 
operation of inverting the chemical source Jacobian matrix can then be 
executed simultaneously for a large number of points on the same vector plane, 
resulting in a large saving of CPU time. 


GEOMETRY AND BOUNDARY CONDITIONS 

A unique subsonic inlet condition procedure for chemically reacting flows 
has been developed in this study for the convergent-divergent nozzle calcula- 
tions. As total mass flow through the thruster is determined by conditions at 
the thruster-throat, the inlet velocitites could not be specified. The total 
enthalpy and mass fractions are specified at the inlet and velocity is extrap- 
olated from the interior. The pressure is obtained using an isentropic 
relationship. 

The subsonic inflow boundary condition procedure is briefly outlined in 
the fol lowi ng . 


Step 1: 

The isentropic assumption at the inflow gives 


dh - 


dp 
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(17) 


which can be written as 
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By assuming frozen species composition, the gas constant, R, is a 
constant, and thus equation (18) can be integrated to yield a relationship 
between temperature, T, and pressure, p, at the inflow, 


T 

CpdT = R In P ( 19 ) 

^ref ^ p ref 


In equation (19) the polynominal function of T for Cp is ultilized 
to obtain an explicit form relating T and p. 

Step 2: 

Step v = 0, and extrapolate u from interior. With known total enthalpy, 
the specific enthalpy can then be calculated. 


Step 3: 

Using equations (4) and (19), temperature and pressure are obtained, and 
density is then calculated from the equation of state. 

For the test problems considered here, the outflows are always supersonic, 
hence, the dependent variables are extrapolated from the interior. Along the 
solid wall, no-slip conditions are specified. The wall is assumed adiabatic. 
Along a plane of symmetry, the normal derivative of all dependent variables 
are zero, except for v velocity, where v = 0 is specified. 

The configuration used for the small thruster simulation is the downsized 
rocket engine design with a large expansion nozzle. A schematic is shown in 
figure 1. Oxygen and hydrogen are used for propellants. A surplus of hydrogen 
was used for wall cooling as well as greater thrust. The thruster inlet has 
two inlet zones. The core flow consists of gaseous hydrogen and oxygen at an 
equivalence ratio of 0.5. The second inlet flow is pure hydrogen used for 
direct wall cooling. Approximately, 5 percent of the outer inlet area is used 
for the second inlet flow. The two inlet flows are separated by a cylindrical 
sleeve with milled slots for the hydrogen. The sleeve and combustion within 
the sleeve are not modeled in this calculation, instead, the calculation starts 
at the point past the sleeve. In the actual thruster combustion is initiated 
by an igniter. The igniter is only used momentarily as the reaction is self 
sustaining. The inlet is rather warm as the propellants were initially used 
for wall cooling. This causes rather high temperatures for the combustion 
products in the core zone. In the calculation the core flow is modeled by a 
3000 < mixture of oxygen and combustion products (table I). The inlet condi- 
tions came from a sequence of initial runs. The first calculation used equi- 
librium mixture. This mixture then changed due to the use of a finite 
rate chemistry model in the code. Mixture fraction values along the 
forward centerline were used for successive calculations. A 1.5 in. long 
combustion chamber is provided where the oxygen rich core flow and wall cooling 
hydrogen can react. The flow is then accelerated through a converging- 
diverging nozzle. The expansion ratio used in the actual thruster was around 
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100. For the present two-dimensional rectangular calculation an area expansion 
ratio of 4.2 is used. The grid for the calculation is shown in figure 2. The 
grid is finer near the walls for improved calculation of viscous effects and 
the mixing and chemical reactions of the shear layer. 

A second calculation is also presented. This involves a large area-ratio 
expansion nozzle. The area ratio is approximately 20:1 and is generated by a 
partial sine curve. The inlet Mach number of 1.1 is specified. The tempera- 
ture and species concentrations are similar to the core flow for the thruster. 
The species concentrations were generated using the equilibium code of Gordon 
and McBride (ref. 12). An inert fraction of 80 percent nitrogen was specified. 
The inlet species fractions are listed in table I. 


RESULTS AND DISCUSSION 
Small thruster 

A contour plot of Mach number for the small thruster calculation is shown 
in figure 3. Mach number increases along fluid flowpath. Mach number contours 
are remarkably similar to other calculations of this thruster geometry, whether 
they involve reaction or not. The maximum Mach number is 2.699 and occurs 
along the exit centerline. Approximately 2/3 of the way through the converging 
part of the nozzle the maximum Mach number (at a fixed axial distance) changes 
from occurring along the centerline to approximately 85 percent of the trans- 
verse distance to the thruster wall. This increase in Mach number at this 
location is caused by the reduction in flow area especially for the outer flow. 

A contour plot of the pressure field is shown in figure 4. Pressure dec- 
reases along fluid flowpath. The highest pressure is found in the area of 
hydrogen injection. Pressure increases slightly with increasing radial 
distance in the converging portion of the thruster nozzle. Pressure is lower 
closer to the thruster wall in the initial expansion. Pressure then becomes 
slightly higher near the walls where velocity is lowest. The lowest pressure 
occurs at the exit centerline (the point of maximum Mach number). 

Contour plots of O 2 , H 2 * anc * H 2 O are shown in figures 5 to 7, 
respectively. As expected, the oxygen mass fraction decreases towards the 
thruster wall, due to the lack of mixing in this region. Oxygen concentration 
does not vary much along the nozzle centerline. The hydrogen mass fraction 
increases with decreasing distance from the thruster wall. The maximum H 2 
mass fraction at the exit is 0.97. The maximum centerline exit O 2 mass 
fraction is 0.455. Much of the oxygen and hydrogen did not mix, causing a 
loss in thruster performance. A turbulence model tailored for a reacting 
shear layer would be more appropriate than the model used for this calculation 
(better mixing would also be obtained by modeling the turbulence generated by 
the core flow combustion and the actual inlet geometry). The water mass 
fraction reaches a transverse maximum where the shear layer is reacting. The 
maximum water fraction is 0.918 and occurs at the exit of the thruster in 
figure 7 (the plotting package truncated the list of contour levels). 

A contour plot of temperature is shown in figure 8. Overall temperature 
drops through the thruster due to expansion and mixing. The maximum radial 
temperature at the inlet shows up at 0.9 of the total radial distance. The 
location of maximum temperature at a fixed axial distance moves further away 
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from the nozzle wall through the combustion chamber and the converging portion 
of the nozzle. The maximum radial temperature then occurs along the thruster 
centerline for the first half of the divergent section of the nozzle. In the 
final half of the divergent section, centerline thermal energy is more 
completely changed to kinetic energy and the location of maximum temperature 
moves away from the centerline. The nozzle wall temperature is much lower 
than the temperature in the mixing layer and the core flow, indicating that 
the hydrogen wall jet has well served as cooling protection for the wall. 

OH mass fraction contours are shown in figure 9. The maximum OH mass 
fraction occurs near the nozzle wall. This does not occur exactly where the 
maximum (sometimes local) radial temperature does, but slightly closer to the 
thruster wall where the hydrogen mass fraction is higher. The location of 
maximum water concentration is also slightly closer to the thruster wall. 


DIVERGING NOZZLE 

The pressure contours are shown in figure 10. Pressure decreases with 
axial distance. The pressure contours are parabolic in shape. This same 
shape is shown by the Mach contours in figure 11, except near the wall, where 
viscous effects predominate. The maximum Mach number for the diverging nozzle 
is 4.103. 

Temperature contours are shown in figure 12. The location of maximum 
temperature is at the nozzle wall, due to viscous heating. In the core flow 
region, temperture initially rises slightly. This is attributed to rec- 
ombination reactions, which release thermal energy to the flow. It is noted 
that the initial large amount of recombination is partly due to the mismatch 
of the equilibrium composition specified at the inflow and the finite-rate 
chemistry model . 

In figure 13, H 2 O contours show a slight decrease at the inlet. This 
is due to the change of temperature and pressure caused by the flow expansion, 
and the change in chemical reaction models, upsetting the inlet equilibrium 
composition. Further downstream from the inlet, the chemical reactions slow 
down and gradually reach the frozen chemistry condition, as evidenced by the 
contour plots. The slowing down of chemical reactions is caused by the 
decrease in temperature and pressure due to expansion. Most of the change in 
water mass fraction with radial distance occurs near the nozzle wall, due to 
the large temperature change near the wall. The hydrogen, oxygen and hydroxol 
concentrations shown in figures 14 to 16, follow temperature contours very 
closely. Oxygen mass fraction initially increases along axial direction, due 
to recombination reactions. Hydroxol mass fraction initially drops to 
50 percent of inlet values and then further decreases with decreasing tempera- 
ture. Oxygen mass fraction increases with decreasing temperature, except for 
a small initial change. Monatomic oxygen mass fraction, shown in figure 17, 
increases rapidly near the inlet showing a mismatch in the equilibrium inlet 
condition and the finite-rate chemistry model in the code. Near the nozzle 
wall, monatomic oxygen mass fraction increases, due to the dissociation reac- 
tions caused by higher temperatures in this region. Monatomic oxygen does 
decrease along the nozzle centerline after the adjustment in chemistry models, 
due to recombination reactions. 
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CONCLUDING REMARKS 


The flowfields in converging-diverging ducts and supersonic nozzles have 
been analyzed using a comprehensive reacting flow code. The calculations show 
the hydrogen wall cooling to be effective for an initial hydrogen-rich thruster 
design. A calculation for a diverging nozzle also shows the code to be capable 
of calculating reacting nozzle flow where the change in area ratio is large. 

The results indicate that the qual i ti tati ve features of the flowfield and 
chemical reactions are well represented by the calculations. However, mixing 
between the wall cooling H 2 jet and the hot gas in the core region for the 
small thruster is underpredicted. A turbulence model more sophisticated than 
the Bladwin-Lomax model used in the present study may be required to correctly 
predict the mixing behavior of the flow. 
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TABLE I. - THRUSTER AND DIVERGING NOZZLE INLET MASS FRACTIONS 
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CONTOUR LEVELS 
0.00000 
0. 10000 
0.20000 
0.30000 
0.40000 
0.50000 
0.60000 

0 . 70000 
0.80000 
0.90000 
1 .00000 

1 . 10000 
1 . 20000 

1 . 30000 

1 .40000 

1 .50000 
1 .60000 
l . 70000 
1 .80000 
1 .90000 
2.00000 
2.10000 
2.20000 

2.30000 

2.40000 

2.50000 
2.60000 
2.70000 



FIGURE 3. - MACH NUMBER CONTOURS. 
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0.65000 
0.70300 
0.75000 
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0.85000 
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0.95000 
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1 . 10000 
1 . 15000 
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LEVELS 



FIGURE 4. - PRESSURE CONTOURS. 


CONTOUR LEVELS 
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FIGURE 5. - OXYGEN MASS FRACTION CONTOURS. 
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CONTOUR LEVELS 
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FIGURE 6. - HYDROGEN MASS FRACTION CONTOURS. 
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0.70000 
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0 . 77500 
0.80000 
0.82500 



FIGURE 7. - WATER MASS FRACTION CONTOURS 
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FIGURE 8. - TEMPERATURE CONTOURS. 
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FIGURE 9. - HYDOXOL MASS FRACTION CONTOURS 




CONTOUR LEVELS 
0. 00000 
0.05000 
0 . 10000 
0. I 5000 
0 . 20000 
0.25000 
0. 30000 
0.35000 
0.40000 
0.4JOOO 
0 . 50000 
0.55000 
0.60000 
0.65000 
0.70000 
0.75000 
0.80000 
0.85000 
0. 90000 
0.95000 
1 . 00000 
I ■ 05000 
1 . 10000 
1 - 15000 



il 


FIGURE ». - PRESSURE CONTOURS. 
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FIGURE 12 


TEMPERATURE contours. 
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FIGURE 13. - WATER NASS FRACTION CONTOURS. 
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FIGURE 14. - HYDROGEN NASS FRACTION CONTOURS 
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FIGURE 15. - OXYGRN MASS FRACTION CONTOURS 
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FIGURE 16. - HYDROXOL MASS FRACTION CONTOURS. 
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FIGURE 17. - MONATONIC OXYGEN MASS FRACTION CONTOURS 
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